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Abstract 

This report presents the formulation and check out problems for 
a computer code DYNAPLAS, which analyzes the large deflection elastic^ 
plastic dynamic response of stiffened shells of revolution. The for- 
mulation for special discretization is by the finite element method 
with finite differences being used for the evaluation of the pseudo 
forces due to material and geometric nonlinearities. Time integration 
is by the Houbolt method. The stiffeners may be due to concentrated 
or distributed eccentric rings and spring supports at arbitrary angles 
around the circumference of the elements. Check out problems include 
the comparison of solutions from DYNAPLAS with experimental and other 
computer solutions for rings, conical and cylindrical shells and a 
curved panel. A hypothetical submarine including stiffeners and 
missile tube is studied under a combination of hydrostatic and dy- 
namically applied asymmetrical pressure loadings. 
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INTRODUCTION 

Due primarily to advances in computer technology the analysis of the 

large deflection elastic-plastic response of realistic structures is now 

within the realm of reality. Several computer codes have already been 

developed for this purpose. These include finite difference codes by 

the groups at M.I.T.,^’^ Ballistic Research Laboratories Sandia 
8 32 57 61 

Laboratories, ’ and Lockheed. ’ Finite element computer codes have 

been developed by McNamara and Marcal , Wu and Witmer, and Stricklin 
49 

et. al. In general all the codes have certain limitations. The 
finite difference codes are generally restricted to unstiffened shells 
whereas the finite element codes do not have this restriction but are restricted 
to simple structures. There is a need for a code to analyze the large de- 
flection elastic-plastic deformation of stiffened shells of revolution. 

Basically there are three different formulations which have been 
used in the analysis of the nonlinear behavior of structures. The first 
formulation treats the effects of nonlinearities as pseudo forces on the 
right hand side of the equations of equilibrium. This formulation re- 
quires that only pseudo forces be computed but has the disadvantage of 
the occurrence of numerical stability problems when certain solution pro- 
cedures are employed. The basic equations of equilibrium neglecting 
damping, for this method are derived in detail in this report and symbol- 
ically may be written as 

[M]{q} + [K]{q} = {P} + {0^} + 


(1) 
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where 

[M] = mass matrix 

[K] = stiffness matrix 

{q} = generalized displacement 

{P} = generalized forces due to applied loads 

{Q^} = pseudo forces due to initial (plastic) strains 
NL 

{Q } = pseudo forces due to geometric nonlinearities 


1 r oc 

Farhoomand and Wilson and McNamara and Marcal use the incremental 
form of Eq. 1 for their formulation. This is obtained as follows, first 
Eq. 1 is written in incremental form: 


[M]{Aq} + [K]{Aq} = {aP} + {AQ^} + {aQ^*-} 


( 2 ) 


I NL 

next the increments of Q and Q are treated as differentials and eval- 
uated by applying the chain rule in terms of the generalized displace- 
ments q^ . 


I ^^i 


,i 


aqNL . . 


qqj = - rNL . 


(3) 


Using the relations given by Eq. 3, Eq. 2 becomes 

[M]{Aq> + ([K] + [K^] + [K'^'-]){Aq} = {aP} + {f^} (4) 

where the unbalance in force {f^} has been added to the right hand side 
as was done in Refs. 18 and 36. 
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A new and completely different formulation is presented by Wu and 
Witmer. Starting with the virtual work expression in terms of stress 
and increments of strain, Wu and Witmer obtained the equilibrium equa- 
tions in the form: 

[M]{q> + {P} + [H]{q> = {F} (5) 

In Eq. 5, {F} represents the generalized forces due to external loads. The 
matrices {P} and [H] depend on the state of stress in the body. 

It is difficult to state the relative advantages and disadvantages of 
the formulations given by Eqs. 1, 4, and 5, respectively. However for the 
asymmetric deflection of shells of revolution where the displacements 
in the circumferential direction are represented by Fourier series the 
representation given by Eq. 1 is superior. This is due to the coupling 
between Fourier terms which appears in [K^], and [H] in Eqs. 4 and 

5. This coupling, for all practical purposes, eliminates the possibility 
of using implicit solution procedures for the formulation given by Eqs. 4 
and 5. 

Regarding solution procedures it should be noted that there is no solu- 
tion procedure which may be designated as "The Solution Procedure" due to 
the dependence of solution procedures on the problem under consideration. 

The objective here is to discuss several solution procedures in general 
and devote special emphasis to the formulation given by Eq. 1 as applied 
to shells of revolution. 

The central difference solution procedure for the time response has 
long been the favorite of researchers using a finite difference formulation 
of the spatial derivatives. 
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More recently, central differences have been used in conjunction with the 

Op 

finite element method by Wu and Witmer, Key and Beisinger, and Krieg 
33 

and Key. The general consensus reached by these researchers is that 
the central difference solution procedure should be used in conjunction 
with the "lumping" of masses at the nodes. Further, Key and Beisinger 
have presented a rational method for lumping the rotary inertia. 

Undoubtedly the central difference or some other explicit solution 
procedure becomes quite attractive as the band width of the stiffness 
matrix becomes reasonably large. It should be pointed out, however, that 
the use of conditionally stable procedures is somewhat a contradiction of 
the basic philosophy of the finite element method. One of the advantages 
of the finite element approach is that the size of the elements may vary 
and small elements may be used in regions where the stress gradient is 
large. However, the time increment which may be used is determined by 
the highest represented frequency of the system which in turn is increased 
by using very small elements. Thus, the time increment for numerical sta- 
bility may become so small as to be of little practical value when very 
small elements are used. This was indeed found to be the case in Ref. 49 
for shells of revolution. 

There are three implicit solution procedures which have received con- 

21 42 

siderable attention. They are the Houbolt, Newmark Beta method, and 
Wilson^ solution procedures. They are similar in that the matrix which 
must be "inverted" is a combination of the mass and stiffness matrix. 

For linear problems the Houbolt procedure is unconditionally stable where- 
as the Newmark Beta method and Wilson procedures are stable for a certain 
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range of parameters. 

6 43 

Bathe and Wilson and Nickell have presented interesting studies 
of the three methods and have shown that the artificial damping is equiva- 
lent to conducting an analysis through modal superposition with the higher 
modes being suppressed. They also present some interesting figures show- 
ing what time increment must be chosen to prevent excessive damping or 
phase shifts. Nickell also presents a discussion of the solution of non- 
linear problems but no numerical examples are presented. 

The authors' experiences with the Houbolt and Newmark Beta procedures 
are reported in Ref. 54 but are worthy of a summary herein. A particular 
form of the Newmark Beta method, commonly referred to as the method of 
Chan, Cox, and Benfield,^^ was used to solve many problems in Ref. 54. 

It was found that both procedures are no longer unconditionally stable 
when geometric nonlinearities are included. However, the Houbolt method 
was found to be considerably more stable than the method of Chan, Cox, and 
Benfield. In both methods the pseudo forces on the right hand side were 
determined based on a linear extrapolation; but it was found that the 
method of Chan, Cox, and Benfield becomes unstable even if the pseudo 
forces are used as their values at the previous time increment. Another 
general conclusion which may be reached from the results presented in 
Ref. 54 is that it is extremely risky to draw conclusions about a non- 
linear analysis based on a study of the linear problem. The same reason- 
ing works in reverse as the method of Chan, Cox, and Benfield is superior 
to the Houbolt method for linear problems. 
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In the present research the Houbolt solution procedure is used al- 
though it is planned to permit an option of either the Houbolt or Newmark 
Beta Methods in the future. The Houbolt method was selected for economic 
reasons but results presented herein show that the method is quite accurate. 
The artificial damping in the Houbolt method permits the use of the econom- 
ical formulation given by Eq. 1 and further permits the pseudo force terms 
to be extrapolated without appreciable loss in accuracy. 


FORMULATION 


There are two basic formulations which have been used in nonlinear 
analysis by the finite element method. The first is the Lagrangian formu- 
lation and the other is the use of an incremental moving coordinate system. 
The three formulations discussed in the introduction which includes the one 
used herein are of the Lagrangian type. The basic starting point for the 
formulation is the equations of equilibrium written in terms of the Kirchhoff 
stress^^ (Fig. 1). 


8U, 


IsT [Sjk(^ik '■o ''oi ' '’o^i 


where 


’jk 


"i 


u- = 


-> 0^0 

«ik 


Kirchhoff stress tensor 
coordinate in original body 
Lagrangian displacement 
body force 
Kronecker delta 


( 6 ) 


u = acceleration 
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6W* = virtual work done in deformed body 
Vq = volume of undeformed body 

The derivation of Eq. 7 from Eq. 6 follows exactly the same procedure 

as for the small deflection case presented by ArgyrisJ Obtaining the 

virtual work in the deformed body, 6W*, requires a physical interpreta- 

45 

tion of the Kirchhoff stress tensor as presented in Novozhilov and, 

1 R 

in more detail, by Haisler. 

Restricting attention to small strains, the Kirchhoff stress is 
the true stress and is related to the elastic component of the Green 
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strains through the matrix [D] 

{S} = [D]{e®} (9) 

As small strains are assumed the total strain is the linear super- 
position of the various components. 

{e} = {e^} + {e^} + {e^} + . . . (10) 


Solving Eq. 10 for the elastic strain, substituting into Eq. 9, and sub- 
stituting the result into Eq. 7 yields: 


/ p u. 
^0 1 


6u^. dV + 


/ (Lej - 






...)[D]{6e} = 6W* 


( 11 ) 


For some problems the potential due to external forces may be a higher 
order function of the displacements; but, as usual, the assumption of 
a first order function of the displacements is assumed herein. 


Thus 


W* = LqJ{P(t)} 


6W* = LPJ{<5q} 


Taking the variation with respect to generalized coordinate q^. 
yields the equation of equilibrium: 

‘ Oo“j 5 " * ^ lIItJ [D]1.) dv 
% ' ''o 


! l|f-J Mf/ldV - / L|f-J [0]{e'^)dV - . 


■aq 


9q. 


( 12 ) 

(13) 


(14) 

{P} 
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It is convenient to write the total strain as 


e = 


=L * 'NL 


( 15 ) 


where and are the linear and nonlinear contributions, respectively, 
to the total strain. Substituting Eq. 15 into Eq. 14 and expanding the 
second term on the left hand side, yields 


^ ^ ^ <^NL> ■*'' 

V % 


/ L 


3e 


NL 


8q. 


•J [D]{£L>dV 


/ Lf^JCDKe'’} dV 


(16) 


/ L|^J [D] 
3q. 


{e'} dV - 


= Pi 


The first term in Eq. 16 produces the terms of the mass matrix times the 
accelerations. The second term gives the contribution to the usual 
linear stiffness matrix times the generalized coordinate. The remain- 
ing terms may be combined to yield 


™ij ‘'j * ''ij ‘’j ' - 


3e 


NL 


/ L^J [D] 


- / l|f7J [D] - •••1 dV (17) 



10 


where 




f P 




0 3q. 


aq. 


dV 


(18) 


"ij ' ^ "'' 

V ‘ ^ 

0 

It should be noted that the volume integrals extend over the entire re- 
gion affected by q. and q.. This integration is, of course, performed 

* J 

by integrating over each element separately and assembling the results 
in the standard manner. 

Writing Eq. 17 for each and every degree of freedom yields the com- 
plete set of equilibrium equations in matrix form. 

[M]{q> + [K]{q} = {P} + {Q*} (20) 

where 

^ [D]{£l} dV - ; [0]{e^L - ...} dV (21) 


The last term on the right side of Eq. 20 is generally called the pseudo 
force and is a function of the unknown displacements. 

In Eq. 21 the pseudo forces due to material and geometric nonlinearities 
are included together instead of separating them into components. The separ- 
ate components are given as: 


Q* = Q? + 0^*- 


( 22 ) 
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where 

= pseudo force due to initial (plastic) strains 
= pseudo force due to geometric nonlinearities 


Q? = / L|^J [D] {e^ + . . .} dv (23) 

1 dq • 

''o 

of = - / L^J [D] (=,} dv - / LffrJ [0] (=nl> 

% ' ''o 

(24) 

9£| 

= - / L^j [D] {=nl) 'i'' - ^ L^j m (=> 

If 1 II 1 


The last form of Eq. 24 is the more efficient from the computational point 
of view when only geometric nonlinearities are considered. Furthermore, within 
the realm of shell analysis Eq. 24 may be integrated through the thickness 
of the shell. The approximate expressions, assuming moderate rotations, for 
total strain in the shell with the ring stiffener being a special case is 
given by: 


{e} = {e} + {e^^} + z{k} (25) 

where 

{e} are the usual expressions for the linear membrane strains 
{k} are the changes in curvature and twist 
and z is the distance from the reference surface. 
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Substituting Eq. 25 into the second form of Eq. 24 and integrating 
through the thickness yields 




i IfrJ 


r L 


ae 


NL 


gq_ -J [D] {e + + zk} dA 


(26) 


where t is the thickness and z is the distance from the reference surface 
to the centroid of the area under consideration. In the present research 
it is assumed that the midsurface of the shell is the reference surface 
and thus z for the shell is zero. But, the circumferential stiffeners 
may be eccentrically located which gives a non-zero value for i. 

The basic governing equations are Eqs. 20, 23, and 26 and these equa- 
tions should be discussed. First from Eq. 23 it should be observed that 
the pseudo forces due to initial strain are functions of the displacement 
and hence vary with time is a second order function of the displace- 
ments). Next it should be observed that the formulation to this point is 
valid for any type of shell. Specialization of the formulation to the shell 
of revolution and ring stiffeners is presented in a later section. 


PLASTICITY RELATIONS 


The Von Mises yield condition and isotropic hardening are used in 
the present study. However, this section includes a discussion of kine- 
matic hardening as well as isotropic hardening. The presentation for 
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38 

isotropic hardening follows that given by Marcal. All discussions are 
for the plane stress case. 

Isotropic Hardening 

The expression for the equivalent uniaxial stress is given by 

a = lal + al - (27) 

where s and e are the meridional and circumferential directions respectively. 
Elastic behavior is observed if a is less than the yield stress in uniaxial 
tension. The normality condition for the increment of plastic strain is: 

{de'’} = df {|^} (28) 

where 

p 

{de } = increment of plastic strain 

_p 

de = increment of equivalent uniaxial plastic strain. 

The hardening rule for the material is simply the relation between 
the uniaxial plastic strain increment de^ and the uniaxial stress increment. 

d5 = H' di^ (29) 


For any type of stress strain curve 


where 

Ey = tangent modulus 
E = Young's modulus 


( 30 ) 
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The other relation needed to complete the formulation is the equation re- 
lating the increment of stress to the increment of elastic strain. 

{da} = [D]{de®} = [D]({dE> - {de^}) (31) 

Premultiplying Eq. 31 by Lt^J and using Eqs. 28 and 29 yields; 

O0 

L||J (da) = d5 = H' de'' = Lffj [D] ({de} - df {|^}) (32) 

_p 

Solving Eq. 32 for de yields: 



H' . L|ij [D] (|i, 


(33) 


In the computational procedure the equivalent uniaxial strain given by 

Eq. 33 is first computed and then the increments of plastic strain and stress 

are computed through Eqs. 28 and 31 respectively. When the equivalent uni- 

-P 

axial strain computed by Eq. 33 is less than zero unloading occurs and de 
is set equal to zero. 

The treatment of the transition range from elastic to plastic behavior 

oo CA 

follows the same procedure given by Krieg and Duffey and Yamada et.al. 

It was found during the course of preliminary research in this area 
that the straight-forward computational procedure presented here can, if 
large increments of strain occur, yield stresses which deviate appreciably, 
from the assumed stress-strain curve. To avoid this deviation the strain 
increment in Eq. 33 is divided into m sub-increments and the procedure 
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repeated for each sub- increment. The stresses are updated in each sub- 
increment which gives modified values for the direction normal to the 
yield surface L|“J* The number of sub-increments is determined by com- 
puting an equivalent uniaxial strain increment, using the relation for 
equivalent uniaxial plastic strain, and dividing by an allowable incre- 
ment of strain. Thus the number of sub- increments varies with time and 
position on the structure. 

The storage requirements for the implementation of isotropic harden- 
ing are the three plastic strain arrays, the uniaxial plastic strain array, 
and the maximum stress array. The uniaxial plastic strain is needed to 
determine H' whereas the maximum stress must be exceeded after unloading 
before additional plastic straining occurs. 

Kinematic Hardening 

The form of kinematic hardening presented here is based on Ziegler's 
modification of Prager's hardening rule.^^ Much of the derivation of 
kinematic hardening follows the same procedure used for isotropic harden- 
ing; but is presented here for completeness. 

The yield condition for kinematic hardening is given by: 

5 = C(o, - ^ (o, - (34) 

where the a's represent the translation of the yield surface. Yielding 
occurs whenever a is greater than the yield stress in uniaxial tension. 

The flow rule for kinematic hardening is, as in isotropic hardening, 
the normality condition. 
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{de^} = df {|^> (35) 

The hardening rule is the one given by Prager 

Ida - H'de^J {|^} = 0 (36) 

where H' is again given by Eq. 30. 

The stress increment is again determined by the elastic strain incre- 
ment through the matrix [D], 

{da} = [D]{de®} = [D]({de} - (de'"} (37) 

Substituting Eq. 35 into Eq. 36 and taking the transpose of the result- 
ing equation yields 


Lfs-J H' dEP 


(38) 


Also using Eq. 35 in Eq. 37 and premultiplying by yields: 

oO 


Lfij (da) = L||j [D]((d.) - d-/ 4)) 
Solving Eqs. 38 and 39 for de*^ yields: 


(39) 


.-P # tde} 

cle = 

L|fJ(C H' J + [D]){||) 


( 40 ) 


The diagonal matrix [ H' J has been inserted for convenience in computa- 
tion. The similarity between Eq. 40 and Eq. 33 should be observed. Using 
kinematic hardening Eqs. 40, 35, and 37 form the basic equations; however. 
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these equations must be supplemented by an equation which yields the incre- 
ments in a^, a-, and a_„. This equation is the primary modification made 
by Ziegler and is given by: 

{da} = dy {a - a} (41 ) 

The term dy which must be greater than zero is determined by noting 
that the yield surface translates but does not enlarge. Thus: 


do = L|fj{da} + L|f-J{da} = 0 
Examination of Eq. 34 reveals that 


(42) 


{— } = 




(43) 


Substituting Eq. 43 into Eq. 42 and using Eq. 41 yields the desired ex- 
pression for dy. 


dy 


L|fj {. - o 


(44) 


Since the stress increment is known from an earlier calculation it is 
a simple matter to evaluate dy through Eq. 44 and then {da} through Eq. 41. 
Thus the computational procedure is complete. 

Kinematic hardening requires seven arrays, three strains, three a's, 
and the equivalent uniaxial strain for the determination of H'. This is 
two more arrays than needed for isotropic hardening. There is no appreci- 
able additional computational effort required to use kinematic hardening. 
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Comparison of Isotropic and Kinematic Hardening* 

Comparing the fundamental equations used for isotropic and kinematic 

hardening reveals several differences. For initial loading the expres- 
-P 

sions for de given by Eqs. 33 and 40 are different. For isotropic harden- 
ing the denominator contains H' whereas for kinematic hardening the same 
term is [|^] [ H' ] It should also be noted that a is defined 

differently in isotropic and kinematic hardening. 

Both isotropic and kinematic hardening have been coded and a compari- 
son made of results. One such comparison for the center deflection of an 
impulsively loaded plate is shown in Fig. 2. It is observed in Fig. 2 
that the results for the initial peak are almost identical for isotropic 
and kinematic hardening. However, the "snap-back" deflection is appreci- 
ably different with kinematic hardening showing that the plate returns to 
the original position. It should be observed that there is appreciable 
strain hardening in the results shown in Fig. 2. To explain this behavior 
consider the simple one dimensional structure shown in Fig. 3, which is 
loaded into the plastic range with appreciable strain hardening and then 
released. Isotropic hardening predicts a final deflection which is approx- 
imately equal to the deflection when the structure is released. However, 
due to the large amount of strain hardening, kinematic hardening predicts 
a plastic instability in the reverse direction. This is obviously a con- 
tradiction to experimental observation and could lead to misleading results. 
It should be noted in Fig. 3 that the value of H' is assumed to depend on 


*Appreciation is expressed to Dr. Nicholas Perrone of the Office of Naval 
Research for his discussion on this subject. 
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the equivalent uniaxial plastic strain only. This could be easily 
changed for a one dimensional state of stress but to the authors' know- 
ledge has not been developed for the plane stress case. 

The net result of this comparison is that neither isotropic nor 
kinematic hardening represents the true situation when appreciable strain 
hardening occurs. However, it is felt that isotropic hardening will yield 
more realistic results for appreciable strain hardening and is used in the 
computer code. 


SOLUTION PROCEDURE 

Equation 20 is solved by the Houbolt solution procedure. The purpose 
of this section is to review the Houbolt procedure and to discuss the approx- 
imation of the loads vector. 

In the Houbolt method the accelerations in Eq. 20 are replaced by a 
finite difference approximation of the second derivative. This substitu- 
tion allows development of recursion relations which can be used for the 
step-by-step calculation of the displacements of the shell and stiffeners. 

The accelerations of the nodes of the shell are approximated by the 
third order backwards difference expression 

- 5% * ''"n-l - %-2> 

2 

The accuracy of this representation is of the order (At) . 

Substituting Eq. 45 into Eq. 20 and simplifying yields 

(2[M] + (At)^ [K]){qn+i> = (At)^ + Q*^^} + [M]{5q^ - 4q^_^ + q^_2> 
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This equation is valid for all time increments but must be modified 
for the first step since the values of {q_ 2 l and {q_-| } are not known. 

This occurrence is common when solving initial value problems by finite 
difference procedures and merely requires that a method be developed to 
start the solution. Equation 46 will, however, be used to calculate the 
displacements for all time increments except the first. 

To start the solution, assume 

“ m 3% - 6q_i + q_ 2 > 

and 

{q.} = {qi - 2q. + q (48) 

° (At)^ ^ ° ^ 

The initial accelerations are obtained using Eq. 20 evaluated at t = 0 
which gives 

[M]{q^} = {P} + {Q*} - [K]{q^} (49) 

Rearranging Eq. 48 gives 

{q_^} = (At)^ {q^} + {2q^ - q^} (50) 

By combining Eqs. 47 and 50 an expression for {q_ 2 > is developed 

{q_ 2 > = 6(At)^ {q^} + 6At {q^} + 9 {q^} - 8 {q^ } (51) 

Substituting Eqs. 50 and 51 into Eq. 46 for the first time increment 

(n = 0) and approximating the forces at time, t = 0, provides an expression 
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in terms of the initial displacements, velocities, and accelerations which 
can be solved to obtain {q-j} 

(6[M] + (At)^ [K]){q^> = (At)^{P^Qj + Q*(o,q)} + [M]{2(At)^ q^+ 6At qg+eq^jl 

(52) 

This equation is used to determine the displacements at the end of the 
first time step. Using Eq. 50, a fictitious matrix of displacements, 

{q_^} can be determined. The displacements are then available so Eq. 46 
can be used for each subsequent time step. 

Examination of Eq. 46 shows that the applied and pseudo loads should 
be evaluated at time increment n+1 in the determination of the displacements 
at increment n+1. However {Q*} is a function of the displacements and thus 
Eq. 46 is basically a nonlinear set of algebraic equations. In the present 
research and as presented in Ref. 54 the problem is linearized by using an 
extrapolation procedure. Further, provisions are made for updating the 
pseudo force vector every m time increments to save on computer time. It 
has been found that if the pseudo loads vector is updated every time incre- 
ment then a full linear extrapolation should be used. However, it has also 
been found that for some problems updating the pseudo loads vector, {Q*}, 
every five time increments and using a half order extrapolation yields ac- 
curate results. For problems exhibiting a high degree of geometric non- 
linearities it is advisable to update the pseudo loads vector every time 
increment to avoid numerical instabilities. 

The option of updating the pseudo loads vector infrequently is quite 
valuable when using the Houbolt solution procedure. If high frequency 
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response is desired it is necessary to use a very small time increment due 
to inherent damping in the basic Houbolt procedure. However, this need 
not cost appreciably in computer run time provided the pseudo forces are 
not updated every time increment. 

The loads vector {P} is determined based on input loads at arbitrary 
times. For intermediate times the loads vector {P} is determined by linear 
interpolation. 


SHELL OF REVOLUTION 

The previous sections have been devoted to structures or shells in 
general. The purpose of this section is to specialize the equations to 
the shell of revolution and present the fundamental theory used in this 
research. 

For the linear stiffness matrix the shell of revolution is idealized 
as a sequence of curved elements. One such curved element is shown in Fig. 

4. The curvature of the element is represented by a second order polynomial 
in the meridional distance s. 

2 

= a^ + a-|S + a2S (53) 

where s is the distance along the meridian of the element. The constants 
aQ, a^ , and a 2 are determined by requiring the actual shell and the idealized 
shell to have the same coordinates and slopes at the nodes. Detailed expres- 
sions for the a's are given in Ref. 53. 
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Strain Displacement Rel at ions 

The strain displacement relations used in this research are based on 

45 

the theory presented by Novozhilov. 

^ 1 
"s = ®s ^ ^^s ^ 2 ®13 

^se " ®se ^'"s0 ^ ®13®23 

where e^, e„, and e^. are the total strains in the meridional and circum- 
ferential directions and the shear strain respectively. It should be noted 
that Eqs. 54 are valid for moderate rotations only. The expressions for 
the various terms in Eqs. 54 are given by: 

A 

= (8U/3S) - ((I'w 

= (l/r)[(8v/30) + u sincfi + w cos(|)] 

o 

~ (l/r)(3u/30) - (v/r) sin(j) + 3v/3s (55) 

So 

e^2 " (3 w/3s) + U((»' 

®23 ~ (3W/30) - (v cos<|))/r 

Xe = - 0/r) (se^j/se) - (1/r) sin* 

A ^ ■ A ^ 

Xs0 = " (1/r) (3e^2/9®) + (sin(|)/r) 


( 56 ) 
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Linear Stiffness Matrix 

The linear stiffness matrix is obtained based on orthotropic material 

and neglecting transverse shear deformations. The strain energy expres- 
sion is given in Ref. 49. The displacement functions used in obtaining 
the linear stiffness matrix are given by: 

W = (a^ + a^S + “ 2 ^^ COS i0 

u = (a^ + OgS + B^s (s-L) + (s-L)) cos ie (57) 

V = (ay + agS + (s-L) + B^s^ (s-L)) sin ie 

Note that the summation convention is being used. L is the meridional 

length of the element. It should be noted that only the terms due to 
symmetric loads about 0 = 0 are included. The terms B-| , ^ 2 * ^ 4 ’ 

are eliminated by static condensation in the Fourier terms i = 0 and 1 
only. These are the only terms in the Fourier expansion which contribute 
to rigid body motion and these terms are needed for that purpose. Static 
condensation is not used for the higher Fourier terms as experience has 
shown that a "too flexible" stiffness matrix may be obtained. 

After the assumed displacement functions and strain displacement 
relations are substituted into the strain energy expression, the element 
stiffness matrix for each Fourier term is determined by numerical inte- 
gration. Twenty-nine Simpson stations are used to integrate along the 
meridian of the shell element. 
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Mass Matrix 

A consistent mass matrix based on linear displacement functions in 
u and V and a cubic displacement function in w is used in this research. 

Rotary inertia is included but has not been found to be significant for 
problems considered to date. 

Nonlinear Terms 

In contrast to the highly accurate procedure used in evaluating the linear 
stiffness matrix, an extremely simplified procedure is used in evaluating 
the effects of both material and geometric nonlinearities. This consists 
of using conical frustum elements and finite difference expressions for 
the strains, rotations, changes in curvature, and twist and evaluating 
the integrals over the length by strip integration. This section presents 
the details of this procedure. 

Referring to the expressions for the strains, rotations, curvatures, 
and twist given by Eqs. 55 and 56 and using the fact that the displacements 
are represented as Fourier series in the circumferential angle e the various 
terms may be written as: 

^ i ^ i 

e_ = e^ cos le e. = e. cos ie 

5 S o o 


" i 

sin 

ie 

" i 


e = p 

se se 



1 - 

sin 


" i 

(58) 

®23 ®23 

ie 

K^ = Kg COS 16 



Kg = k] cos ie = <1^ sin ie 
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where e^, e^, e^g, ej^, e 23 > k^, k^, and <^g may at most depend only on 
the meridional distance, s, along the element. However, using finite 
difference expressions for these terms their values may be computed at 
the middle of the element and held constant over the meridional length. 

The use of finite difference expressions permits the individual 
Fourier components to be written as: 





J 


de 


aq. 


se 

r^j 


ae 


aq 


kqi 

1 


'13 


ae 


aq 


6 


ae 


aq 


3k 


aq 




(59) 


1 

^0 = 




9q. 


9k 


'se 


9q 




There is no summation on i, which indicates the particular term in the 
Fourier expansion. The detailed expressions for the partial derivatives 
are given in Appendix A. The partial derivatives in Eq. 59 are independent 
of s, e, and the displacements, but may depend on the Fourier number. 
Specializing Eq. 23 to the shell of revolution yields: 


= fff L"^^j[D]{e^ + + ...} rdodsdz 


(60) 



where i and j are the Fourier number and degree of freedom respectively. 
All terms are assumed to be constant over s, then 


/ ds = L 

Further, the expression is simplified with the following definition 

{C} = [D]{e^ + J + ...} rL 


where r = radial distance to the middle of the element. 

Using Eqs. 54, 58, and 59 the expression for the pseudo force may be 
written as 


Ji 


9e^ . aeL 3 k^ 

ff [ — ^ cos ie + e,- — ^ cos ie + z — cos ie) C, 


3e' . 9eL 3 k' 

+ — J- COS ie + e«- — ^ sin ie + z — f cos ie) C, 

3,1 233,1 3,1 


ae' 


ae. 


ae. 


”se ^ °''i '3 

+ ( — ^ sin ie + e-j^ — sin ie + 623 — ^ cos ie 


aq. 


aq. 


aq. 


9k 


+ z — sin ie) C~] dz de 
aq. 


where there is no summation on i. C.|, C 2 , and are the three components 
of {C} given by Eq. 62. 
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Applying the same logic to Eq. 26 the pseudo forces due to geometric 
nonlinearities for the shell may be written as: 

Ml •; 

Q. = - / [-4 cos ie C, + — ^ cos ie C. + sin ie C.] de 
3q. aqj aq. 


9e]- i . ael- 

f [e-|2 — 4 ^1 ®23 — ^2 


aq. 


aqj 


(64) 


* (e 




9e«« ^ 

sin ie + e 


^ cos ie) C-] de 

9q 


where there is no summation on i and 


{C} = (t)(r)(L) [D] {e^L> 


(65) 


{C} = (t)(r){L) [D] {e + 

Equations 63 and 64 are the basic equations needed to compute the 
pseudo forces due to material and geometric nonlinearities. The C's de- 
pend on z and e in Eq. 63 and C and C depend on e in Eq. 64. The integral 
through the thickness is performed by Simpson's rule in Eq. 63 and it has 
been found that 5 to 7 Simpson stations are quite adequate. The integrations 
around the circumference in Eqs. 63 and 64 are performed using a modified 
Simpson's rule, the modification being that the second order function of e 
is weighted with sin ie and/or cos ie to obtain the integrals. Before 
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this modified Simpson's rule was incorporated it was found that the number 
of stations around the circumference was dictated by the number needed to 
integrate cos ie or sin ie and could be quite large. 

RING STIFFENERS 

The ring stiffeners for any element may be constructed of as many 
as three different segments and may vary from element to element. This 
permits an exact representation of stiffeners in the form of T or I 
sections. However a Z section may be represented only approximately as 
the product of inertia terms are neglected. 

The mass, stiffness, and strain hardening for the stiffeners may be 
different from those of the adjoining shell. The effects of eccentricity are 
accounted for in all calculations where the reference surface is taken 
to be the mid-surface of the shell. The ring stiffeners are assumed to 
be in a state of uniaxial stress in the circumferential direction. Under 
this assumption strain energy per element based on linear theory is given 
by 

L EN ^ ^ p ^p 

Ur = ^ / [eg + Z e0 K0 + Kg] rdedzds (67) 

where N is the number of stiffeners per element. Integrating Eq. 67 
through the thickness yields 

yL = / [^2 + j gg <0 + I Kg] rdeds (68) 

where z is the distance from the mid-surface of the shell to the centroid 
of the stiffeners and I is the area moment of inertia about the mid-surface 
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of the shell . 

The terms of the stiffness matrix are obtained from Eq. 68 by taking 
the second derivative with respect to the generalized degrees of freedom. 


l<R - L 
ij 3q^-9qj 


(69) 


In the numerical computation of the contribution of the stiffeners to the 
element stiffness matrix the N stiffeners are assumed to be "smeared" over 
the meridional length of the element. Integration is performed along the 
length using Simpson's rule. 

For most practical applications the ring stiffeners are discrete and 
"smearing" can lead to appreciable error. This error is avoided by using 
a very short element for the discrete stiffeners, thus avoiding any "smear- 
ing" error. 

The pseudo forces due to initial (plastic) strains are given by Eg. 

63 with only being non zero. Thus the pseudo force expression becomes: 




N ff [( 


3e 

3q^ 


cos 10 + Z 




3q. 


COS i 9 ) E ( + 


. . ) rL] dz de 

(70) 


Equation 70 is evaluated using either strip or Simpson integration over 
each segment of the stiffener and a modified Simpson intearation around 
the circumference. Note that strip integration has been used over the 
meridional length. 

The pseudo forces due to geometric nonlinearities are given by Eq. 
26 which when specialized to the ring becomes: 
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qNLi 


-2 

- (E)(t)(r)(L)(N) f [( — ^ cos ie + z — ^ cos ie) %eL] de 

aql aq] 


( 71 ) 


ae 


- / e 


23 


'' aqV 


sin ie (e^ + % 023 ^^ + z k^) de 


where z is the distance from the mid-surface of the shell to the centroid 
of the complete ring stiffener. There is no summation on i in either 
Eqs. 70 or 71. 


SPRING SUPPORTS 

There are frequently shells of revolution which have intersecting 
supports at various locations around the circumference. The missile 
tubes and platforms in submarines serve as specific examples of such 
supports. In the present research these supports are included as linear 
elastic springs and are incorporated into the basic equations as pseudo 
forces on the right hand side of the equations of equilibrium. This pro- 
cedure is rather straightforward but can lead to numerical instability 
problems for overly stiff supports. Further the supports must be in- 
cluded so as not to prevent rigid body translation of the complete 
system. To accomplish this latter objective the supports are treated as 
coupled linear elastic springs. 

The general equations for the forces and moments due to the spring 
supports are given by: 

{Fg} = [Kp] {q} (72) 

3x1 3x33x1 
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and 

{M} = [Km] (73) 

where q are the deflections at the center of the element in cylindrical 
coordinates, , ^ 2 , and q^ are the axial circumferential and radial dis- 
placements respectively. e ^3 and e 23 are rotations about the circumfer- 
ential and meridional directions respectively. 

As a simple example consider a single support on a circumferential 
ring as shown in Fig. 5. The support as shown acts in the q-j direction 
only and the forces are given by: 



where q^ and q^ are in the other two orthogonal directions. 

Next the q's are related to the q's through the transformation 



Using Eq. 75 to transform the displacements and forces F', the matrix 
[Kp] becomes 
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[Kp] = 


0 0 0 


cos0g 0 0 


sin0g 0 0 


K 0 0 


0 0 0 


0 0 0 


0 COS0 sin0 


0 0 


0 0 


Carrying out the indicated matrix multiplication in Eq. 76 yields: 


(76) 


[Kp] = K 


0 


0 cos 0, 


0 cose^sine^ 
. s s 


cos0^sin0^ 
s s 


sin^e. 


(77) 


The matrix of coefficients given by Eq. 77 restricts motion in the q^ ' 
direction but permits the entire system to translate in the horizontal 
direction. It should be noted that this objective could not be accom- 
plished without permitting a complete matrix of spring constants. 

An interesting application of the spring support idea is to repre- 
sent an incomplete clamped circular ring or cylindrical panel. The spring 
constants in these cases should be very large - - ideally infinity. How- 
ever in practice these spring constants may be selected by assuming that 
the supports are short stubby beams as compared with the actual structure. 
It has been found that the spring constants may be made quite large with- 
out introducing numerical stability problems. Specific examples are given 
in the section on application. 
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COMPUTATIONAL PROCEDURE 

The purpose of this section is to present, in some detail, the com- 
putational procedure used to compute the pseudo forces due to geometric 
and material nonlinearities. This should not be interpreted to imply 
that the other portions of the numerical process are unimportant. How- 
ever, the details of much of the procedure for computing mass and stiff- 

49 53 

ness matrices are the same as were presented in earlier works. * 

In an earlier computer code called DYNASOR for the geometric non- 
linear analysis of shells of revolution the integrals around the 
circumference were evaluated exactly. To perform this task several three 
and four dimensional arrays were needed where each dimension was the num- 
ber of Fourier terms used in the expansion. Further, these arrays were 
used in three and four nested DO loops within the program. The net result 
of the formulation used in DYNASOR was that storage requirements and com- 
puter run times restricted the number of Fourier terms that could be used. 

During the course of the present research, it was found that numerical 
integration around the circumference is a very efficient computational pro- 
cedure both from the viewpoint of computer time required and storage al- 
locations. The importance of the modified Simpson's rule cannot be over- 
emphasized. The present code called DYNAPLAS is considerably more effici- 
ent especially when a large number of Fourier terms is being used. Con- 
sidering the fact that DYNASOR required shorter computer run times than 
comparable computer codes, the present code, DYNAPLAS, should be capable 
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of solving large scale problems in realistic computer run times. 

The pseudo forces are computed in one rather long (1400 statement) 

subroutine called QPRIME. The -long length is due to the large number 
of statements inside a DO loop over the number of modified Simpson sta- 
tions around the circumference. This also explains why a reduction in 
the number of stations around the circumference through the use of a 
modified Simpson's integration saves appreciably on computer run times. 

The key to an efficient computational procedure is the computation 
of the partial derivatives given in Eq. 59 and tabulated in Appendix A. 

While the term partial derivatives is used here they could equally well be 
called coefficients or any similar terminology as they do not depend 
on the displacements. At times the term partial derivatives leads one 
to believe that the Eqs. 59 are only approximate whereas they are exact 
regardless of the magnitude of the generalized displacements. The second 
step in the computational procedure is to compute linear strains, rotations, 
changes in curvature, and twist for each Fourier term as given by Eq. 59. 
With these preliminary calculations out of the way a DO loop over the num- 
ber of integration stations around the circumference is entered. For each 
modified Simpson station the following calculations are performed. 

A. Compute the strains, rotations, curvature, and twist in accordance 
with Eq. 58. 

B. If stations have not yielded before, check for yielding. This cal- 
culation involves calculating the stress in accordance with Eqs. 9 
and 10 and the evaluation of a from Eq. 27. 

C. For each Simpson station through the thickness compute the increment 
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of plastic strain. This involves Eqs. 28> 30, 31, and 33. Further, 
this section involves some rather complex logic to transverse the 
region from elastic to elastic-plastic behavior and checks for un- 
loading and reloading. 

D. Compute the pseudo forces due to initial strains using Eg. 63. 

E. Compute the pseudo forces due to geometric nonlinearities using Eg. 64. 



COMPUTER PROGRAMS 


The computational algorithm presented here has been coded into two 
FORTRAN programs called SAMMSOR III and DYNAPLAS. These two codes are 
essentially an extension of two SOR codes, SAMMSOR II and DYNASOR II, 
which have been operational since 1970. The extensions include elastic- 
plastic material behavior, ring stiffeners, and the effects of other 
internal and/or external stiffening members in addition to the large 
deflection capability of DYNASOR II. 

As in the SOR series of codes the dynamic analysis is conducted 
by first executing the SAMMSOR III code to obtain an output tape con- 
taining the stiffness and mass matrices for the particular geometry 
being studied. DYNAPLAS is then executed to solve the dynamic equa- 
tions for a specific set of initial conditions, boundary conditions 
and loading history. A subsequent analysis of the same problem with, 
for example, a different loading history requires only the execution 
of DYNAPLAS. 

The SAMMSOR III code utilizes a highly refined curved shell of 
revolution element in addition to beam type ring stiffeners. The shell 
element utilizes cubic displacement functions for the normal and in- 
plane displacements and, through static condensation, a basic eight 
degree of freedom stiffness matrix is generated. A mesh generating 
routine is included which allows the user to input the geometry of 
complicated shells of revolution with only a minimum of input informa- 
tion. The ring stiffeners are assumed to have zero products of inertia 
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but may be eccentrically located and may be formed by as many as three 
rectangular flange members. 

In addition to the stiffness and mass matrices generated by 
SAMMSOR III, DYNAPLAS requires the uniaxial stress-strain data for the 
shell and stiffeners, the boundary conditions, the initial displacement 
and velocity conditions, and the applied load as a function of time. 

In addition, a number of other control constants are required; for ex- 
ample the specific Fourier terms to be utilized in the analysis, the 
time increment, the number of integration stations to be utilized in 
numerically integrating over the volume, print parameters, etc. The 
applied load history may be described by specifying either a pressure 
distribution over the element or the Fourier coefficients of the distri- 
bution at discrete time intervals with a linear variation being assumed 
between the specified times. For a particular element, the pressure 
distribution is assumed to be constant in the meridional direction and 
vary as a step function in the circumferential direction. In addition, 
the code is capable of accepting concentrated ring loads at each node. 
The uniaxial stress-strain data is described by inputting a piecewise 
linear curve with as many as five stress-strain data points. The ma- 
terial characteristics may vary from element to element. The code may 
of course be run with plasticity and/or geometric nonlinear effects 
omi tted. 

The computer code has a restart provision permitting the program to 
be restarted at a particular time increment once the program has been 
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run up to that time. Restart information is placed on tape periodically 
during the execution of DYNAPLAS. If subsequent analysis of the output 
by the user indicates that, for example, a smaller time increment is 
needed or the nonlinear forces need to be updated more often, the pro- 
gram can be automatically cycled to any time increment for which restart 
information is stored on tape and then the analysis restarted with new 
input parameters (smaller time increment, etc). 

The program output consists of all input control words, input loads, 
generalized forces, stiffness and mass matrices, and the resultant dis- 
placements, stresses and strains through the thickness and around the 
circumference. Various print parameters allow the user to select as 
much or as little of the above output as he desires. 

In order to make the program more flexible variable or dynamic 
dimensioning has been used on forty-four main arrays in DYNAPLAS. A 
blank common block is dimensioned for a certain number of storage loca- 
tions in the main program and is used as a dynamic storage area. The 
forty-four arrays are then implicitly equivalenced to various portions 
of the blank common by the subroutine call list. These variables are 
variable dimensioned in all subroutines for the number of elements, 
harmonics, etc. A subroutine is called from the main program which scans 
the input files to determine the necessary storage requirements based 
on the number of elements, harmonics, etc. If the required number of 
storage locations is less than what the blank common is dimensioned, then 
execution of the data set begins. This subroutine also determines (based 
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on the input data) the beginning point of each of the forty-four arrays 
relative to the first location of the blank common. This feature allows 
the user to set the size of the blank common based on (a) available core 
size and (b) program size required to solve a problem. The values of 
the various input parameters (such as number of elements, Fourier terms, 
integration stations, etc.) may be input in any combination as long as 
the dynamic storage area is not exceeded. 

These programs have been written in standard FORTRAN IV language. 

The programs are operational on the IBM 360/65 computer at Texas A&M 
University and the CDC 6600 computer at Sandia Laboratories. As usual, 
the IBM version of the programs requires double precision arithmetic 
in critical areas. Utilizing mixed mode arithmetic, SAMMSOR III re- 
quires approximately 180,000 bytes of storage (fifty elements). DYNAPLAS 
requires approximately 160,000 bytes of storage plus the dynamic storage 
area. A fifty element idealization utilizing ten Fourier terms requires 
approximately 200,000 bytes of dynamic storage area (IBM 360/65). 

Considering the complexity of the computer program it is highly 
efficient. For example the solution to the ring presented in Figure 2 
required less than 20 seconds of IBM 360/65 computer time. The analysis 
of the circular plate presented in Figure 9 required about 3 minutes of 
IBM 360/65 time. The solution to the asymmetrically loaded truncated 
cone shown in Figure 11 using thirty elements and ten Fourier terms 
required 30 minutes of CDC 6600 computer time (updating the nonlinear 
forces at each time increment). If the nonlinear forces are updated 
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every three time increments, the solution may be obtained in only 10 
minutes of CDC 6600 computer time. 

The SAMMSOR III code has been operational for several years and 
its validity demonstrated on many test cases. The DYNAPLAS code has 
been in operation for about one year at Texas A&M University and 
Sandia Laboratories. Based on the test cases reported herein, good 
agreement has been noted between DYNAPLAS and numerical results ob- 
tained using other computer codes. Considering the assumptions in the 
plasticity theory, reasonable agreement has been achieved between 
DYNAPLAS and experimental results. 

Copies of the users manuals* and computer code are available to 
qualified users upon request. Requests should be addressed to Dr. Walter 
E. Haisler, Aerospace Engineering Dept., Texas A&M University, College 
Station, Texas 77843. 


*1 . Haisler, W.E., Stricklin, J.A., and Von Riesemann, W.A., "DYNAPLAS - 
A Finite Element Program for the Dynamic, Elastic-Plastic, Large 
Deflection Analysis of Stiffened Shells of Revolution," SLA-73-0127, 
Sandia Laboratories, Albuquerque, New Mexico (Also Rpt. 72-27, Aero- 
space Engineering Department, Texas A&M University), January 1973. 

2. Haisler, W. E. , Stricklin, J.A., and Von Riesemann, W.A., "SAMMSOR 
III - A Finite Element Program to Determine Stiffness and Mass 
Matrices of Ring-Stiffened Shells of Revolution", SLA-73-0126, 

Sandia Laboratories, Albuquerque, New Mexico (Also Rpt. 72-26, 
Aerospace Engineering Department, Texas A&M University, College 
Station, Texas), January 1973. 
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Check Out Problems 

The purpose of this section is to present the solutions to a sub- 
stantial number of problems. The problems vary from a simple symmetri- 
cally deformed ring to the response of a hypothetical but complex sub- 
marine. The solutions show the versatility of DYNAPLAS as well as its 
limitations. 

Static Solution for Spherical Cap 

A doubly curved shell element is used in the computation of the 

linear stiffness matrix, but, conical frustum elements are used in 

the calculation of the effects due to nonlinearities including the 

53 

stress resultants. As past experience has shown that the use of 
only conical frustum elements gives a residual bending moment, it was 
deemed necessary to check the accuracy of using conical frustum ele- 
ments for the computation of the stress resultants. 

The problem chosen to check the accuracy of the computation of 
the stress resultants was the static solution for a linearly elastic 
spherical cap under an internal pressure as shown in Fig. 6. The cap 
has a radius of 100", a thickness of .5", a Young's modulus of 10 x 10^ 
psi, and a Poisson's ratio of 0.2. The static solution was obtained by 
using a large time increment in the solution procedure. The large time 
increment introduces a considerable amount of artificial damping in the 
basic Houbolt solution procedure and thus yields a static solution 
after a relatively few time steps. This same procedure is used on 
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the submarine problem, presented later, to obtain the solution under 
a hydrostatic pressure before the transient load is applied. 

The results shown in Fig. 6 are for TO and 30 equally spaced ele- 
ments in DYNAPLAS and the converged solution. It is noted that the 30 
element idealization is in excellent agreement with the converged solu- 
tion. Further there is no residual bending moment for either the 
coarse or fine element breakdown. The conclusion drawn from this 
study is that the stress resultants may be accurately computed based 
on conical frustum elements. 

Finally it should be explained why conical frustum elements are 
used for nonlinear terms. The answer is simply that the use of 
finite difference expressions to compute the strains using curved 
elements shows straining under rigid body motion. 

Symmetrically Impulsively Loaded Ring 

The first example of an elastic-plastic dynamic response is a 

plane strain ring subjected to a syrranetrical impulse loading giving 

an initial velocity of 4911.7 in/sec. The ring had a 10" radius, a 

thickness of .1", and a density of .1 Ib/in? The yield stress was 

taken to be 42,400 psi and the problem was studied for various amounts 

of strain hardening. The results for the radial deflection vs. time 

are shown in Fig. 7 along with results given by Duffey and Krieg^^ 
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and computer results obtained from the computer code UNIVALVE. 

Results are presented for three different values of the strain 
hardening parameter, x = Ej/E. It is noted that all three solutions 



are in excellent agreement. A problem similar to this one was solved 
where the impulsive load was applied as a high intensity pressure over 
a short duration of time. The results were found to agree with those 
obtained using an initial velocity based on same impulse. 

Free Ring Under Half- Cosine Impulsive Load 

The first example solution for asymmetrical loading is a free- 
free plane strain ring under a half-cosine impulsive loading. The 
aluminum ring had a radius of 11", a thickness of .15", and a density 

3 

of .09997 Ib/in . The stress-strain curve was assumed to have a yield 
stress of 30,000 psi, a secondary modulus of 5 x 10^ psi to a total strain 
of .009 in/in, and perfectly plastic thereafter. The solution through DYNA- 
PLAS was obtained using 5 Fourier terms and a time increment of 0.5 y 
sec. The results for the deflection at e=0 vs. time are presented 

O 

in Fig. 8 along with results from the computer codes SPECTRE and 
3 ? 

UNIVALVE. Results from DYNAPLAS are given for updating the pseudo 
forces every time increment in conjunction with an extrapolation 
factor of 1.0 and updating every 5 time increments with an extrapola- 
tion factor of .5. It is noted that the two solutions agree quite 
well with each other as well as with the results from SPECTRE and UNIVALVE. 

This problem was used to check various provisions in the computer 
code. The problem was solved as a ring, a shell (with slight modifi- 
cation to computer code) and as three flanges, two of which were 
eccentrically located. The problem was also solved as a ring with 
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the reference surface at the inner edge of the ring. The solutions 
for all these cases agreed quite well with each other. 

Figure 9 presents the plastic strain at the outer surface and 
0 = 0 vs. time for 3, 5, and 9 modified Simpson's stations around 
the circumference. It is noted that reasonably good results are 
obtained with only 3 modified Simpson's stations. The deflections 
agreed quite well with the results presented in Fig. 8 for all three 
cases. 

Clamped Ring Under Impulsive Loading 

Wu and Witmer ’ have conducted extensive studies of the highly 

nonlinear dynamic response of rings and reported experimental results 

based on earlier tests. One of the rings which they studied is shown 

in Fig. 10. The ring had a radius of 2.935 in., a thickness of .123 

-3 2 4 

in, and a density of .25 x 10 lb-sec /in . The yield stress was 
42,800 psi and the secondary modulus was 78,700 psi. The ring sub- 
tended an angle of 315° and the boundaries were clamped. An impulse 
giving an initial velocity of 4862 in/sec was applied over a 120° 
segment. 

The theoretical solution by DYNAPLAS was obtained using 5, 10, 
and 15 Fourier terms. Time increments used were 5 y sec for the 5 
and 10 term analysis and 2 y sec for the 15 term analysis. It should 
be noted that the deflections are very large with the outer edge of 
the ring almost reaching the origin of the original circle. The 15 
Fourier term solution shows excellent agreement with the experimental 
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results. 

The clamped boundary condition was obtained by including springs 
at the supports. In an earlier run, one of the spring constants was 
arbitrarily set equal to 15 x 10^ Ib/in. It was found that this very 
large spring constant introduced a "flip-flop" response in the solu- 
tion near the support. The portion of the ring very near the spring 
would be loading plastically in one direction at one time increment 
and be loading plastically in the opposite direction at the very next 
time increment. As this provision is unrealistic and not accounted 
for in DYNAPLAS, the stress deviated from the assumed stress-strain 
curve. However, the solution procedure did not diverge. To avoid 
this difficulty a physically reasonable support system was assumed 
to be represented by the bending of a strip with the same cross- 
section as the ring but only .5" long. This led to spring constants 
of 50,000 Ib/in for the deflection springs and 25,000 in-lb/rad. for 
the rotational spring. The deflections and rotations at the boundary 
were so small as to be insignificant for all computer runs using 
these spring constants. The important point to be noted is that 
physically realistic spring constants must be used. 

The case for 10 Fourier terms was first studied using 11 modi- 
fied Simpson stations to integrate around the circumference. This 
led to a monotonically divergent solution at approximately 1000 y 
sec. The problem was corrected by increasing the number of modified 
Simpson stations to 17. This is the only problem where this type of 
difficulty occurred but as a result it was concluded that a "safe" 
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number of Fourier stations is approximately twice the highest 
Fourier number being used. Undoubtly most problems may be solved 
using considerably fewer stations and thus saving appreciably on 
computer run times. 

In the solution using 15 Fourier terms, 27 modified Simpson 
stations around the circumference and 7 Simpson stations were used 
through the thickness. The pseudo forces were updated every time 
increment resulting in a computer run time of 8 minutes (CDC 6600). 

The results from DYNAPLAS for the total strain at the line of 
symmetry and on the outer surface are given in Fig. 11. 

It is noted that the maximum strain reaches a value of 8.5%. 

The angle of rotation as a function of the angle around the circum- 
ference at 1000 y sec is shown in Fig, 12. It is noted that the 
maximum rotation is .6 radians at 1000 u sec and probably approaches 
1 radian at 2000 y sec. The rotations were obtained from the solu- 
tion using 10 Fourier terms. 

Cylindrical Shell Under Impulsive Load 

The final deflection of a clamped cylindrical shell under a 
half-cosine impulsive loading is given in Fig, 13. The length of 
the cylinder was 6 inches and was idealized using 10 finite elements 

along half the length and 5 Fourier terms. The final deflection was 
obtained by running DYNAPLAS with a small time increment until all 
elements began unloading. The time increment was then increased 
and the solution restarted. This large time increment damped out the 
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motion quite rapidly. Shown in Fig. 13 are results from DYNAPLAS, 

57 35 

SHORE, and experimental results. Reasonable agreement among 
the results is noted but it should be pointed out that the final 
deflection shape is not a good measurement of the accuracy of a 
computer code. 

Circular Plate Under Impulsive Load 

The next problem is the analysis of a circular plate shown in 
Fig. 14 and 15. The experimental results were taken from Ref. 15. 

In the experimental procedure the plate was supported by two massive 
rings 2" in thickness and connected with 12 1/2 inch bolts at a radius 
of 4 inches. The edge of the rings was 3" from the center. It was 
observed after the experiment that some slippage did occur at the 
inner edge of the ring. To represent the boundary conditions in the 
theoretical analysis it was assumed that the circular olate was 
clamped at the bolt circle and on rollers for the portion inside 
the bolt circle covered by the ring. These boundary conditions are 
duplicated in the sketch in Figs. 14 and 15. 

In the analysis using DYNAPLAS eleven finite elements were used. 
Nine were equally spaced over the inner portion of the ring and two 
elements were used to represent the portion of the ring under the 
rollers. A time increment of 1 p sec was used in DYNAPLAS. The 
stress-strain curve was represented by an initial elastic modulus 
of 10^ psi to a strain of ,00424 in/in, then a secondary modulus of 
3.29 X 10^ psi to a total strain of .00449 in/in, and perfectly 
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plastic behavior thereafter. 

Experimental and theoretical results for the deflection at the 
center of the plate are shown in Fig. 14. It is observed that there 
is good agreement through the initial peak deflection and over part 
of the unloading curve. However the experimental results deviate 
from theoretical results at later time increments. 

The experimental and theoretical values for the meridional 
strain on the top of the plate at a distance of 2.125 inches from 
the center are shown in Fig. 15. Considering the closeness of this 
point to the edge of the supporting ring the agreement between theory 
and experiment is considered to be acceptable. Further if the strains 
were plotted to a scale necessary to represent the strain at the center 
of the plate very little difference in the experimental and theoretical 
points would be observed. 

Truncated Cone Under Half-Cosine Impulse 

The next problem studied was the large deflection elastic-plastic 

dynamic response of a truncated cone under a half cosine impulsively 

applied pressure. The truncated cone had a upper radius of 7.95", 

a lower radius of 10.23", a thickness of .5430" and a density of 
-4 2 4 

1.88 X 10 lb sec /in . Additional details are given on Figs. 16, 

17, and 18. The material was assumed to have a yield stress of 
30,000 psi and to be elastic-perfectly plastic. 

Results are presented for the deflections (Fig. 16) and strains 
(Figs. 17 and 18) as obtained from DYNAPLAS, REPS I and SHORE. 
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Considering the fact that the three computer codes are completely in- 
dependent the agreement of displacements and strains presented in 
Figs. 16, 17, and 18 is considered to be outstanding and serves as 
a check on the accuracy of all three codes. However, based on the 
shown displacements, it is clear that the degree of geometric non- 
linearity is not severe. 

In DYNAPLAS, the conical frustum shell was idealized as 30 
equally spaced finite elements and 10 Fourier terms were used. Seven 
Simpson stations were used through the thickness and a 2 u sec time 
step was used. Two runs were made varying the frequency of updating 
the pseudo forces and the number of modified Simpson stations around 
the circumference. Results were the same for both cases. In the 
first 13 Simpson stations were used and the pseudo forces were up- 
dated every three time increments with an extrapolation factor of 
1.0. The computer run time was 10 minutes on the CDC 6600. In the 
second run 17 Simpson stations were used and the pseudo forces were 
updated every time increment. The computer run time was 30 minutes 
on the CDC 6600. Storage requirements were 88,000 words for the 
second case. 

The SHORE code was run using 30 and 18 equally spaced increments 
along the meridian and around the circumference respectively. The 
computer run time was 22-1/2 minutes on the UNIVAC 1108 which is com- 
parable to the case one run for DYNAPLAS. 
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Cylinder Under Moving Pressure Load 

To illustrate the capabilities of DYNAPLAS to treat a complex 
moving pressure load, a clamped cylinder previously solved through 
REPSIL was chosen. The shell had a length of 5.958", a diameter of 
5.958", a thickness of .042", and a density of .1 Ib/in . Young's 
modulus and Poisson's ratio were 10.5 x 10^ psi and .3 respectively. 
The yield stress was chosen to be 44,000 psi and the material was 
taken to be elastic-perfectly plastic. 

The pressure loading is given by 


p(re,t) 


where 


r 0 for t < ^ 




^ + T-t) cos ^ for 


0 for t > ~ + T 


ire<f 


0 for ~ < 1 re 1 < 


D = diameter = 5.958" 

U = velocity of pressure front = 24,800 in/sec. 

T = duration of triangular pulse = 8.28743 x 10"^ sec. 
Pq = peak pressure at crown line = 28,000 psi. 


The wave front velocity requires 200 p sec to travel to e = tt/ 2. The 
short duration of the pressure pulse yields a very high intensity pressure 
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extending over a very small angle of the shell. 

In DYNAPLAS the shell was idealized as 10 equally spaced ele- 
ments over 1/2 of the meridional length. The Fourier coefficients 
for the pressure loading v/ere input every 4 y sec with the computer 
code using linear interpolation at intermediate times. A time in- 
crement of 1 y sec was used with the pseudo forces being updated 
every time increment. 

The first run was conducted using 5 Simpson stations through 
the thickness and revealed a shortcoming of the DYNAPLAS computer 
code. The rotations at the nodes diverged with the value at one 
node increasing in the positive direction and the adjacent node 
diverging in the negative direction. After considerable deliberation 
it was concluded that the difficulty was caused by the extreme thin- 
ness of the shell (.042"). The elastic stiffness matrix simply has 
no resistance to rapid changes in angles. It was decided that for 
extremely thin shells it should be permissible to use effectively 
membrane theory to evaluate the effects of nonlinearities. This 
consists of using strip integrations through the thickness based 
on one point at the midsurface of the shell. 

Three separate runs were conducted using DYNAPLAS. The first two 
were with 15 Fourier terms, 27 modified Simpson stations around the 
circumference but with and without the rotational degree of restraint 
at the fixed end. If membrane theory is adequate the rotational de- 
gree of restraint should not be important which was in fact found to 



be the case. The 20 Fourier term analysis was conducted as the 
Fourier coefficients for the pressure load were converging very 
slowly. However, it should be recognized that the stiffnesses of 
the higher Fourier terms are quite large and thus a convergent series 
for the pressure loads is not necessary. For 15 Fourier terms and 
27 stations around the circumference the computer run time was 53 
minutes on the IBM 360/65. This is equivalent to approximately 13 
minutes on the CDC 6600. 

The results for the radial deflection and the meridional strain 

at the middle of the shell and at e = 0 are shown in Figs. 19 and 

24 

20 along with the results from REPSIL. Needless to say the agree- 
ment is rather poor. More interesting is that DYNAPLAS predicts 
lower deflections but higher strains than predicted by REPSIL. To 
check the consistency of strains and deflections an elementary 
analysis was conducted. The deflection shape was assumed to be 
represented by 

0 ) - 2 (1 ” cos 

where 6 is the deflection at the center. The average value of the 
strain was then computed by evaluating the deformed length, substract- 
ing it from the undeformed length, and dividing by the undeformed 
length. The results based on the deflections obtained from DYNAPLAS 
are shown in Fig. 20. A check of the strains computed by DYNAPLAS 
reveal that the value of 20 % is very close to the average value at 
the maximum load. Using the value of 2.5" for the maximum deflection 
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as obtained by REPSIL in Eq. 78 yields an average strain of 34%. 

There are two differences in DYNAPLAS and REPSIL which may account 
for the differences in results for this highly nonlinear problem. First 
DYNAPLAS permits moderate rotations whereas REPSIL permits arbitrarily 
large rotations. Second REPSIL takes the pressure as acting normal to 
the deformed surface of the shell where DYNAPLAS uses the reference 
surface of the undeformed shell. This second point may be quite im- 
portant for this problem. Underwood has agreed to study this problem 
and his results should be available before the results are published 
in the open literature. 

Conical Frustum Under Half -Cosine Impulse 

Since it was necessary to use partial membrane theory on the 
previous problem it was deemed necessary to show that the difficulty 
is caused by the thinness of the shell and not because of highly 
nonlinear behavior. For this purpose the previously studied conical 
frustum shell was again analyzed with the initial impulse being in- 
creased by a factor of 2.5. Results for deflections and strains are 
shown in Figs. 21, 22, and 23. Examination of the output data reveals 
deflections over 1/2 of the average radius of the shell and strains as 
great as 40%. Seven Simpson stations were used through the thickness 
and a detailed study of the data revealed no oscillations in the ro- 
tations at the nodes. 
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Hypothetical Submarine 

To illustrate the application of DYNAPLAS to a problem which 
involves a shell, stiffeners and springs, a structure with a vague 
resemblance to a submarine was chosen. Three difference types of 
response are of interest. The first is due to a hydrostatic pressure 
followed by a high intensity transient load. The final stage of re- 
sponse is the overall vibration of the vehicle over a substantial 
period of time. 

Except for the fluid-structure interaction, the total response is 
comparable to what might be encountered by a complex structural system, 
such as a submarine vehicle, subjected to hydrostatic pressures and 
shock environments. The application would prove extremely useful in 
the design and isolation of internal equipment subjected to shock 
loadings of such severity as to cause gross plastic deformation (but 
not complete failure) of the submarine hull. 

In this example, the vehicle is 360" in length, 60" in diameter 
and has an intersecting vertical tube 12" in diameter. The material 
properties are typical of steel and are given in detail in the User's 
Manual for DYNAPLAS. The purpose of the present presentation is to 
present an overall view of the analysis procedure including some typical 
curves. 

Ideally, the analysis would consist of inputting the detailed 
geom.etry of the complete vehicle and conducting the analysis. However 
this would require excessive computer run time and more storage than 
is available on most computers. For this reason the portion of the 
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vehicle around the tube was selected as the region of primary interest. 
The vehicle was idealized as shown in Fig. 24 with the element breakdown 
shown in Fig. 25. Six elements were used for the hemispherical nose and 
rather large elements were used for the stiffened cylinder except around 
the tube where 19 rather closely spaced elements were used as shown in 
Fig. 25b. It was assumed that the vehicle is symmetric about the center- 
line. A total of 45 finite elements were used in the idealization. 
Details of the geometry including the stiffeners are shown in Fig. 24. 

The spring constants for the tube were determined based on elemen- 
tary beam theory. These included both deflection and rotational springs. 
The detailed values for the spring constants are given in the User's 
Manual . 

The hydrostatic solution for a uniform pressure of 500 psi was 
obtained by applying a linearly increasing pressure to a time of 20,000 
ysec and holding it constant at 500 psi to 50,000 ysec. A time step of 
1000 usec was used in the solution procedure. A slow application of 
the pressure was necessary as it was found that a step loading caused 
numerical instabilities when large time steps were used. These in- 
stabilities are believed to be due to the presence of the tube. Also 
it was necessary to apply a radial restraint at the apex to eliminate 
rigid body motion. It was found that the deflection increased linearly 
to 20,000 usec and remained constant thereafter. A check of the results 
shows a membrane solution for the hemispherical nose but variations of 
the stress near the stiffeners and the tube as would be expected. Five 
Fourier terms were used in the solution. 
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Next the boundary condition was removed and a high intensity 
pressure with a decay time of 40 ysec was applied to the shell in 
addition to the hydrostatic pressure. The pressure was assumed to 
vary as a half cosine around the circumference and to decrease as it 
approached the end of the vehicle. For this part of the analysis the 
program was restarted with a time step of 2 ysec and allowed to run 
for 120 ysec. 

For the transient response it was assumed that the only region 
where plasticity and geometric nonlinearities were important was the 
area around the missile tube. A very high yield stress was used for 
the other elements and stiffeners and only a single modified Simpson 
station was used around the circumference. 

To account for the overall response of the vehicle the code was 
restarted at 50,120 ysec with a time increment of 5 ysec and run to 
50,320 ysec. Next it was restarted with a time step of 15 ysec and run 
to 51,400 ysec. In a separate restart at 50,320 ysec a time step of 
50 ysec was used and the analysis was conducted to 55,000 ysec. 

Results for the meridional strain near the tube are shown in 
Fig. 26. The distinct phases of the response to each of the different 
loads can be observed. First the hydrostatic value is reached and then 
a rapid increase occurs when the transient load is applied. The slight 
dip in the value of the strain when the transient load is applied is 
unexplained. The rest of the response is the vibration of the modes of 
the vehicle. 

Two fundamental modes of vibration of the overall vehicle were 
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observed. The first is a lateral vibration similar to a free-free 
beam. This is depicted in Fig. 27 as the response of Fourier coefficient 
No. 1. The actual results showed appreciable lateral translation of the 
vehicle. 

The other response which was easily observed is the response of 
the shell in a ring mode of vibration. An elementary calculation re- 
vealed that this mode has a period of 10,000 ysec. The computations 
were not carried out to where this mode reaches its peak. 

The purpose of this problem was to demonstrate that DYNAPLAS is 
capable of solving such complex problems and not to report a detailed 
analysis due to the assumptions in structure and loads. However it is 
concluded from this study that DYNAPLAS may be used for such analyses. 

Cylindrical Panel 

A cylindrical panel under an impulsively applied load was studied 
to determine if such problems can be solved by representing the bound- 
aries by spring supports. The panel circumscribed an angle of only 
60°. Analyses were conducted using 5, 10, 15, and 20 Fourier terms. 

Results agree with experimental and other theoretical results up 
to near the peak deflection but DYNAPLAS results show considerable more 
"snap-back." A check of results shows that there are appreciable de- 
flections in the shell beyond the spring supports. The conclusion 
reached from this study is that DYNAPLAS may be used for the initial 
response of panels but may not be used to study complete panel behavior 
unless the panel circumscribes almost 360° as was the case for the 
clamped ring under impulsive loading. 



Users Hints 


Recent convergence studies have shown that most results given 
in this report were not obtained in the most efficient manner. The 
following guidelines should save appreciably on computer run time. 

1. For moderately nonlinear problems the psuedo forces should 
be updated every two or three time increments and an extrapolation 
factor of 1.0 should be used. For problems involving a high degree of 
geometric nonlinearities the psuedo forces must be updated every time 
cycle. 

2. The number of modified Simpson's stations around the circum- 
ference may be less than the number of Fourier terms used. 

3. Three to five Simpson's stations through the thickness should 
be adequate for most problems. 



Extensions in Progress 


60 


DYNAPLAS is currently being extended to provide the following; 

1. Variable thickness in the circumferential direction. 

2. Strain rate effects, 

3. Temperature effects including the variations of material 
properties with temperature. 

4. Improved plasticity relations including orthotropic plasticity. 
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Statement of Policy 

The primary objective of this research is to develop a computer 
code which will be useful to the engineering community. Consistent 
with this objective the authors welcome comments from the users and 
are willing to assist the users within reason. This is especially 
true of users associated with the Navy, AEC, and NASA who have sup- 
ported and continue to support this research and the many users of 
SNASOR II and DYNASOR II who have helped, through their comments, 
to make these better computer codes. 
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